close all
clear
%clc
addpath(genpath(pwd))
warning off
%% Section 1: User settings for the estimation
startYear     = 1961;       %Start year for estimation
endYear       = 2019;       %End year for estimation

lambdaSMM     = 1           %Weight in objective function to Shrinkage of based on Method of Moments (SMM)
SMMwithCSreg  = 1           %2 for incl. CS and risk-adj CS moments,
                            %1 for incl. CS moments,
                            %0 for not including CS moments

appOrder      = 2;          % Approximation order of DSGE model
appOrderVar   = [appOrder;appOrder;appOrder]; % Used appOrder per variable: c_t, pi_t, evf_t
endoMeanGrid  = 0;          % 1 for endogenously recenter the grid given the projection approx.
                            % 0 for using the mean in Perturbation to re-center the grid 
PFsolve       = 0;          % 1 for solving under perfect foresight, 0 for normal RE solution

muGrid        = 2;          % mu in Smolyak grid
gridFactor_x  = 1.5;        % grid for proj using +/- gridFactor_x*std1st approx
nx            = 6;          % #of states
ne            = 5;          % #of shocks

numCPUs       = feature('numcores'); % Number of CPUs

numOptimSteps = 0;          % Number of optimizations
optim         = 1;          % 1 Nelder-Mead.
                            % 2 Nelder-Mead with transformed parameter space
                            % 3 Own gradient based optimizer
                            % 4 CMAES
                            % 5 fmincon
cmaesPopSize  = 48;         % population size for CMAES
cmaesSigma    = 0.01;       % sigma in CMAES
cmaesSigmaMax = 0.1;        % sigmaMax in CMAES
MaxIter       = 5000;       % Max number of iterations for the optimizer
MaxEvals      = 10000;      % Max number of function evaluations for the optimizer
TolFun        = 1D-3;       % Function tolerance at the objective function
TolX          = 1D-3;       % Function tolerance for the coefficients
epsValue      = 1e-5;       % stepsize for computing gradient in own optimization routine
seOn          = 1;          % 1 for computing standard errors, else 0
ScalingMacro  = 0.2;        %Setting measurement errors for macro variables
ScalingYield  = 0.1;        %Setting measurement errors for yields

if numCPUs > 8
    % We use Grendel-S
    rmpath([pwd,'/Mex_myPC'])
    MEXon     = 1*0;          %1 For using MEX-files, else 0
else
    % We use My PC
    rmpath([pwd,'\Mex_GrendelS'])
    MEXon     = 1*0;
end

% Default setting
tau           = 400*0;      %Moments computed based on T*tau simulations, else closed-form solution used when tau = 0.
pruningOn     = 1;          %1 for using the pruned approximation, 0 for the unpruned approximation
nyMom         = 10;         %number of the first nyMom observables used for the SMM part
CSregress_m   = 4;          %2 for biannual implementation of CS regressions, 4 for annual etc.
inclSurveys   = 0;          %1 for including surveys, else 0
%% Section 2: Data and Model parameters
% Setting the calibrated parameters as in the sample ending in 2007Q4
structData = loadData_2020(startYear,endYear);
meanData   = nanmean(structData.data,1)';
setCalibrateParams

CSdata     = CampbellSchillerRegression(structData.yieldsAll(1:end,:),CSregress_m);
%[CSdata.maturities;CSdata.CSbetta(2,:); CSdata.CSBetta_CI95Boot]
%[CSdata.maturities;CSdata.CSbetta(2,:);(CSdata.CSbetta(2,:)-1.96*CSdata.CSbetta_se(2,:));(CSdata.CSbetta(2,:)+1.96*CSdata.CSbetta_se(2,:))]
PredData    = predictiveRegression(structData.yieldsAll(1:end,:),CSregress_m);

% Compute empirical moments and weigthing matrix
dataSMM = momentsSMM(structData.data,structData.yieldsAll,nyMom,CSregress_m,SMMwithCSreg);
qLag = 3;
W    = getOptimalWeighting(qLag,dataSMM);

% The variables which we do not select for estimation are calibrated to the assigned values
calibrateParams.CRRA    = 10;

% Starting values
params.ALFA    = -20;
params.BETTA   = 0.995;
params.B       = 0.1;
params.PHI     = 0.075;
params.CHI     = 2.8;
params.XI      = 0.75;
params.KAPAw   = 0;
params.ETA     = 6;
params.PHIpai_1= 2;    %forward looking taylor - rule used
params.PHIc_1  = 0.0;
params.RHOz    = 0.98;
params.RHOd    = 0.99;
params.RHOn    = 0.95;
params.RHOp    = 0.99;
params.RHOa    = 0.99;
params.PAIss   = 1.0;
params.SIGMAmuz= 0.0004;
params.SIGMAd  = 0.03;
params.SIGMAn  = 0.03;
params.SIGMAp  = 0.001;
params.SIGMAa  = 0.01;


% With select = [2:10]*4/CSregress_m in SMM
calibrateParams        = rmfield(calibrateParams,'U0'); %use this to have U0 being active
if lambdaSMM == 0 && SMMwithCSreg == 0
    calibrateParams.ScalingYield = 0.083;
    paramsValues = [ -9.022776063819    0.982049504405    0.842649534371    0.156529406162    5.034468977177  ...
        0.496248691492    0.834187560183    9.989927144903    4.761601472044    4.567610947850  ...
        0.542385704275    0.970250320252    0.982626774885    0.985227451255    0.980480403080  ...
        1.033674162070    0.004028142279    0.151415431980    0.087000694341    0.002472833532  ...
        0.006768681145 ]';
    %Q =  -42.123, done
    %params = rmfield(params,'ETA');
    %calibrateParams.ETA = 10;
    %paramsValues = [  -9.019709489925    0.988082906359    0.843087970113    0.156456240000    5.034964799813  ...
    %    0.496137277258    0.834634654189    4.761693997169    4.566755444927    0.542433116416  ...
    %    0.970196259553    0.983025331472    0.984845599620    0.985596364745    1.034228932383 ...
    %    0.004078909062    0.151410865842    0.087047038405    0.002533510660    0.006769633771  ]';
    % Q = -42.2862, done

    if inclSurveys  == 0
        paramsValues = [ -9.028716590975    0.982266250604    0.842361732428    0.156479081156    5.034836503333  ...
            0.496759474262    0.832590190898    9.993210930796    4.763078785386    4.566576746778  ...
            0.540078967798    0.969923377313    0.983383802732    0.986490376269    0.979508181507  ...
            1.033563585772    0.004029330530    0.152200577559    0.087013513509    0.002469217148  ...
            0.006769472137]';
        %Q =  -33.5685 , done
        params = rmfield(params,'ETA');
        calibrateParams.ETA = 10;
        paramsValues = [  -9.031556169776    0.987786064725    0.840693890552    0.156475393296    5.041790544842  ...
            0.496914550034    0.833584636632    4.763543292616    4.571640015301    0.540083221597  ...
            0.968417548873    0.984048551793    0.985474506211    0.983839479366    1.034411511677 ...
            0.004029158115    0.152130215328    0.087017750871    0.002469224360    0.006773357049 ]';
        % Q =  -34.1646, done

    end
end

if lambdaSMM == 1 && SMMwithCSreg == 0
    calibrateParams.ScalingYield = 0.083;
    paramsValues = [   -9.028038895774    0.998951492708    0.573446962099    0.092142363688    7.109745043309  ...
        0.309835564590    0.907272586310    9.999976637060    3.689998023448    0.657908266214  ...
        0.372438917126    0.985014676053    0.984646051130    0.989483972724    0.981580745194  ...
        1.020839316381    0.003778277331    0.122342761981    0.064388867447    0.000595024634  ...
        0.002181313644 ]';
    % Q =  -35.3757, done
    %params = rmfield(params,'BETTA');
    %calibrateParams.BETTA = 0.999;
    %params = rmfield(params,'ETA');
    %calibrateParams.ETA = 10;
    %paramsValues = [ -9.053582808205    0.576811489694    0.092293743651    7.090051953600    0.309536824749  ...
    %    0.908440605488    3.708774453952    0.657893962864    0.372393424263    0.984945084021  ...
    %    0.984612868642    0.989419283213    0.982013532390    1.020796571634    0.003783169506  ...
    %    0.122558055211    0.063198213207    0.000599056555    0.002202224156 ]';
    % Q = -35.5727, done 
    if inclSurveys  == 0
       paramsValues = [  -9.028029271246    0.998950921316    0.573446350764    0.092142424750    7.109737619856  ...
            0.309835355067    0.907272242017    9.999987557349    3.689999903014    0.657907564838  ...
            0.372438520080    0.985013987989    0.984651325634    0.989483433573    0.981580853441  ...
            1.020839785755    0.003778273304    0.122342737863    0.064388820413    0.000595024423  ...
            0.002181314409 ]';
        % Q =  -27.0467, done
        params = rmfield(params,'BETTA');
        calibrateParams.BETTA = 0.999;
        params = rmfield(params,'ETA');
        calibrateParams.ETA = 10;
        paramsValues = [   -9.028863169865    0.573766013580    0.092007885442    7.104235632347    0.309343776128  ...
            0.909187735984    3.691356617240    0.657924791062    0.372686810520    0.984935427939  ...
            0.984914308540    0.988870881500    0.981422619750    1.020752587707    0.003783245472  ...
            0.122238974768    0.064467237333    0.000596850613    0.002182832852   ]';
        %Q = -27.1532, done 
    end

end

if lambdaSMM  == 1 && SMMwithCSreg  == 1      
    calibrateParams.ScalingYield = 0.083;   
    % Best specification
    paramsValues = [-11.101942823607    0.988637503556    0.635426561457    0.089931047515    5.382449932453  ...
        0.380779218722    0.909884449017    8.801097086628    3.360386068027    1.034982531614  ...
        0.469664278984    0.984459501458    0.989512812916    0.995163270115    0.983275042959  ...
        1.023564607272    0.003720077294    0.128367989135   0.040220146959    0.000407800759  ...
        0.002139052407  ]';
    % Q =   -22.2404, done

    if inclSurveys  == 0
        paramsValues = [-11.460211161892    0.990618453585    0.652915418149    0.097705260203    5.827819204624 ...
            0.437879445269    0.917091392607    8.848256715573    3.957100237398    1.096397434567...
            0.490843444232    0.983587090508    0.986125716788    0.995315021181    0.982541178986...
            1.020706348956    0.003295178790    0.127137250725    0.043709077256    0.000430616002 ...
            0.002363438615 ]';
        % Q = -14.3984, done
    end
end
selectParams = fieldnames(params);
params       = values2struct(paramsValues,selectParams);

%% Section 3: Setting setupFilter and setupProj
% The lower and upper bounds for the parameters
[lowerBounds,upperBounds,Insigma] = setLowerUpperBounds;
upperBounds.BETTA = 0.999;

% The information needed for projection
Smolyak_elem_iso                  = Smolyak_Elem_Isotrop(nx,muGrid);
setupProj.Smol_grid_iso           = Smolyak_Grid(nx,muGrid,Smolyak_elem_iso);
[n_nodes,epsi_nodes,weight_nodes] = Monomials_2(ne,eye(ne));
if PFsolve ==1
    % to get perfect foresight solution
    n_nodes = 1;
    epsi_nodes = zeros(1,ne);
    weight_nodes = 1;
end

% Constructing the struct setupFilter
constructSetupFilter
setupFilter.zlbON          = 1;
%% Section 3: Estimation
tic
[Q0,errorMes,model0,~,outFilter0] = objectFuncQML(struc2values(params,setupFilter.selectParams),setupFilter);
Q0
toc

if numOptimSteps > 0
    for i=1:numOptimSteps
        % We switch between the optimizer in optim and the Nelder Mead
        if mod(i,2) == 1
            setupFilter.optim = optim;
        else
            setupFilter.optim = 6;
        end
        outEstim = runOptimization(setupFilter,params);
        params = outEstim.paramsOpt;
        printToScreen(struct2array(params))
    end
else
    setupFilter.optim = 0;
    outEstim = runOptimization(setupFilter,params);
end

%% Model diagnostics
% The nominal short rate
rAnnual_ss        = outEstim.model.params.PAIss/(outEstim.model.params.BETTA*outEstim.model.params.MUZss^(-outEstim.model.params.CHI*(1-outEstim.model.params.CHI0)-outEstim.model.params.CHI0))^4-1
rAnnual_mean      = outEstim.outFilter.resMoments.meanYields(1,1)
outEstim.outFilter.Q
outDiagnostics = plotOutputFilter(outEstim);
plotTermPremia(outEstim)
plotCSloadings(CSdata,outEstim);
reportFit
CRRA = getCRRA(outEstim.model.params)
outEstim.newsDecomSim = inflVarianceRatio(outEstim,outEstim.outFilter.outSim.xAll_cu);
outEstim.newsDecom    = inflVarianceRatio(outEstim,outEstim.outFilter.xhat);
outEstim.newsDecomSim.summaryTable
outEstim.newsDecom.summaryTable
outEstim.nonMarketUtility = getNonMarketUtility(outEstim);
%error('stop')

% Correlation of the preference shock with the sentiment index by the
% University of Michigan. 
%close all
%figure(10)
%hold on
%plot(zscore(outEstim.setupFilter.sentimentMich),'-k','LineWidth',2)
%plot(zscore(outEstim.outFilter.xhat(3,:)),'-ok')
resSent = nwest(zscore(outEstim.setupFilter.sentimentMich),[ones(length(outEstim.outFilter.xhat(3,:)),1) ...
    zscore(outEstim.outFilter.xhat(3,:)')],4);
[resSent.beta(2) resSent.beta(2)-1.96*resSent.se(2)  resSent.beta(2)+1.96*resSent.se(2)]

%% Additional ex post estimation analysis
outEstim.UnconditionalMoments = unconMoments(outEstim.outFilter.ySimAsInData,outEstim.setupFilter);
outEstim.reducedForm          = doReducedFormAnalysis(outEstim);

if PFsolve ==0
    % The timing premium
    pai0         = 0.05;
    numSimTiming = 2000;
    pathLength   = 5000;
    outEstim.paiOpt = GetTimingPrem(outEstim,pai0,numSimTiming,pathLength);
end

% Computing GIRFs
if seOn == 1
    close all
    IRFlength = 20;
    shockSizeScalar = 1;
    xsGIRF             = [mean(outEstim.outFilter.xhat(7,:)); zeros(5,1)];
    outEstim.GIRFs     = getImpulseResponses(outEstim,IRFlength,shockSizeScalar,[],xsGIRF);

    plotGIRFsChart(outEstim.GIRFs,2)

    xfGIRFxhat         = outEstim.outFilter.xhat(1:6,:);
    xsGIRFxhat         = [outEstim.outFilter.xhat(7,:);zeros(5,setupFilter.T)];
    outEstim.GIRFsxhat = getImpulseResponses(outEstim,IRFlength,shockSizeScalar,xfGIRFxhat,xsGIRFxhat);
    xfGIRFsim          = outEstim.outFilter.outSim.xAll_cu(1:6,1:1000);
    xsGIRFsim          = [outEstim.outFilter.outSim.xAll_cu(7,1:1000);zeros(5,1000)];
    outEstim.GIRFsSim  = getImpulseResponses(outEstim,IRFlength,shockSizeScalar,xfGIRFsim,xsGIRFsim);

    plotGIRFsChart(outEstim.GIRFsxhat,1)
    plotGIRFsChart(outEstim.GIRFsSim,5)
end

testOfSpanningON = 1;
if testOfSpanningON == 1
    numPCAs             = 5;
    outEstim.macroSpanData_5 = doMacroSpanRegressions(structData.yieldsAll(1:end,:),structData.data(1:end,:),setupFilter,numPCAs);
    outEstim.macroSpanData_6 = doMacroSpanRegressions(structData.yieldsAll(1:end,:),structData.data(1:end,:),setupFilter,numPCAs+1);

    % To address comment from referee
    outEstim.macroSpanData_nonLin_5 = doMacroSpanRegressions_nonLinear(structData.yieldsAll(1:end,:),structData.data(1:end,:),setupFilter,numPCAs);
    outEstim.macroSpanData_nonLin_6 = doMacroSpanRegressions_nonLinear(structData.yieldsAll(1:end,:),structData.data(1:end,:),setupFilter,numPCAs+1);
    
    % On Simulated data with Measurement errors
    yieldsAllSimWithNoise        = outEstim.outFilter.yieldsAllSimWithNoise;
    ySimAsInDataWithNoise        = outEstim.outFilter.ySimAsInDataWithNoise;
    outEstim.macroSpanSim_5      = doMacroSpanRegressions(yieldsAllSimWithNoise,ySimAsInDataWithNoise,setupFilter,numPCAs);
    outEstim.macroSpanSim_6      = doMacroSpanRegressions(yieldsAllSimWithNoise,ySimAsInDataWithNoise,setupFilter,numPCAs+1);
    
    outEstim.macroSpanSim_nonLin_5      = doMacroSpanRegressions_nonLinear(yieldsAllSimWithNoise,ySimAsInDataWithNoise,setupFilter,numPCAs);
    outEstim.macroSpanSim_nonLin_6      = doMacroSpanRegressions_nonLinear(yieldsAllSimWithNoise,ySimAsInDataWithNoise,setupFilter,numPCAs+1);

    % On Simulated data withOUT Measurement errors
    yieldsAllSim = outEstim.outFilter.yieldsAllSim;
    ySimAsInData = outEstim.outFilter.ySimAsInData;
    outEstim.macroSpanSimNoNoise_5 = doMacroSpanRegressions(yieldsAllSim,ySimAsInData,setupFilter,numPCAs);
    outEstim.macroSpanSimNoNoise_6 = doMacroSpanRegressions(yieldsAllSim,ySimAsInData,setupFilter,numPCAs+1);

    outEstim.macroSpanSimNoNoise_nonLin_5 = doMacroSpanRegressions_nonLinear(yieldsAllSim,ySimAsInData,setupFilter,numPCAs);
    outEstim.macroSpanSimNoNoise_nonLin_6 = doMacroSpanRegressions_nonLinear(yieldsAllSim,ySimAsInData,setupFilter,numPCAs+1);

    % a linear model - without measurement noise
    modelLin         = outEstim.model;
    modelLin.gxx     = modelLin.gxx*0;
    modelLin.gxxx    = modelLin.gxxx*0;
    modelLin.Gxx     = modelLin.Gxx*0;
    modelLin.Gxxx    = modelLin.Gxxx*0;
    modelLin.hxx     = modelLin.hxx*0;
    modelLin.hxxx    = modelLin.hxxx*0;
    modelLin.Hxx     = modelLin.Hxx*0;
    modelLin.Hxxx    = modelLin.Hxxx*0;
    modelLin.byxx    = modelLin.byxx*0;
    modelLin.Byxx    = modelLin.Byxx*0;
    modelLin.byxxx   = modelLin.byxxx*0;
    modelLin.Byxxx   = modelLin.Byxxx*0;
    if inclSurveys  == 1
        modelLin.Rxx     = modelLin.Rxx*0;
        modelLin.Rxxx    = modelLin.Rxxx*0;
    end
    setupFilterLin   = outEstim.setupFilter;
    setupFilterLin.pruningOn = 0;
    
    % Moments computed by simulation
    setupFilterLin.numSim = 1D5;
    outSimLin      = simProjection(modelLin,setupFilterLin.shocks(:,1:setupFilterLin.numSim));
    if setupFilter.pruningOn == 1
        yieldsSimLin  = getYields(modelLin,outSimLin.xAll_cu);
    else
        yieldsSimLin  = getYields(modelLin,outSimLin.x_cu);
    end
    pos_l             = find(strcmp({'$l_t$'},modelLin.labely));
    pos_w             = find(strcmp({'$w_t$'},modelLin.labely));
    pos_c_ba1         = find(strcmp({'$c_{t-1}$'},modelLin.labelx));
    pos_myz           = find(strcmp({'$\mu _{z,t}$'},modelLin.labelx));
    pos_c             = find(strcmp({'$c_t$'},modelLin.labely));
    pos_pai           = find(strcmp({'$\pi_t$'},modelLin.labely));
    ySimAsInDataLin   = zeros(setupFilterLin.numMacro+length(modelLin.maturities),setupFilterLin.numSim);
    ySimAsInData(1,:) = outSimLin.y_cu(pos_l,:)-modelLin.gSS(pos_l,1);
    ySimAsInData(2,:) = outSimLin.y_cu(pos_w,:)-modelLin.gSS(pos_w,1);
    ySimAsInData(3,:) = 4*(outSimLin.y_cu(pos_c,:)-(outSimLin.x_cu(pos_c_ba1,:)+modelLin.gSS(pos_c_ba1,1))+outSimLin.x_cu(pos_myz,:)+log(modelLin.params.MUZss));
    ySimAsInData(4,:) = 4*outSimLin.y_cu(pos_pai,:);
    ySimAsInData(setupFilter.numMacro+1:setupFilter.numMacro+length(modelLin.maturities),:) ...
                         = 4*yieldsSimLin(modelLin.maturities,:);
    yieldsAllLin        = 4*yieldsSimLin;

    outEstim.macroSpanSimNoNoiseLin_5 = doMacroSpanRegressions(yieldsAllLin,ySimAsInData,setupFilter,numPCAs);
    outEstim.macroSpanSimNoNoiseLin_6 = doMacroSpanRegressions(yieldsAllLin,ySimAsInData,setupFilter,numPCAs+1);
end
% Do table for nonlinear spanning
toTableMacroSpanNonLinear(outEstim)

% Deleting some simulated series
outEstim.outFilter = rmfield(outEstim.outFilter,'yieldsAllSim');
outEstim.outFilter = rmfield(outEstim.outFilter,'yieldsAllSimWithNoise');
outEstim.outFilter = rmfield(outEstim.outFilter,'ySimAsInData');
outEstim.outFilter = rmfield(outEstim.outFilter,'ySimAsInDataWithNoise');


% Standard naming of output file
filename = ['Model_u_endYear',num2str(endYear),'_lambdaSMM',num2str(lambdaSMM),'_SMMwithCSreg',num2str(SMMwithCSreg),'_inclSurveys',num2str(inclSurveys)];
if seOn == 1
    filename = [filename,'_withSE'];
    [struct2array(outEstim.paramsOpt)' struct2array(outEstim.paramsSE_eps5.qLag3.paramsSE)' struct2array(outEstim.paramsSE_eps6.qLag3.paramsSE)' struct2array(outEstim.paramsSE_eps7.qLag3.paramsSE)' ]
end
if params.SIGMAmuz == 0
    filename = [filename,'_SIGMAmuz0'];
end
if params.SIGMAd == 0
    filename = [filename,'_SIGMAd0'];
end
if params.SIGMAn == 0
    filename = [filename,'_SIGMAn0'];
end
if params.SIGMAp == 0
    filename = [filename,'_SIGMAp0'];
end
if params.SIGMAa == 0
    filename = [filename,'_SIGMAa0'];
end
if PFsolve ==1
    filename = [filename,'_PFsolv'];
end
save(filename,'outEstim','outDiagnostics','structData','CSdata')
